1
ECEN 4532: DSP Laboratory
Professor Michael Perkins
Lab 5
March 5th, 2018
2
Contents
1. Introduction 3
2. The JPEG Compression Algorithm 3
3. Conclusion 27
A. MATLAB Code 27
A.1. MATLAB code used to compile report 27
A.2. MATLAB functions 29
A.2.1 dctmgr.m 29
A.2.2 idctmgr.m 30
A.2.3 quantization.m 30
A.2.4 inverseQuantization.m 30
A.2.5 zigZag.m 31
A.2.6 inverseZigZag.m 31
A.2.7 run_bits_value.m 33
A.2.8 inverse_run_bits_value.m 34
3
1 Introduction
In this lab we are implementing the JPEG algorithm used for image compression. We do this so
that images can be transferred and stored efficiently and to save space. Uncompressed images
are very difficult to store, so storing these images digitally is very advantageous.
We are using an encoding/decoding algorithm in order to achieve compression on luminance
(gray-level) images. The main part of this lab is using the Discrete Cosine Transform (DCT) to
describe an image using a smaller set of coefficients to represent the original image. These values
are put through a quantization table to decrease the information needed to represent the
coefficients created through the DCT. Then we decode the information from these steps to get a
new image that is compressed that uses less space than the original picture.
2 The JPEG Compression Algorithm
To implement the JPEG compression algorithm we start by taking the forward and inverse DCT. To
implement the forward DCT the original image is divided into non overlapping blocks size 8x8.
Each block gets transformed using the 2D DCT expressed as,
These coefficients are an approximation and do introduce distortion and slight loss in image quality.
Given the coefficients 𝑓
𝑢,𝑣
̃
one can use the inverse DCT to reconstruct an estimate of the block using
4
1) The function dctmgr is included at the end of the report in section A.2.1.
2) The function idctmgr is included at the end of the report in section A.2.2.
The next step in the JPEG standard is quantization. The effect of the DCT is to create many small
coefficients that are close to zero, and a small number of large coefficients. The goal of
5
quantization is to represent a continuum of values with a finite small set of symbols. In JPEG the
quantization of the DCT coefficients other than the DC coefficient is performed as follows.
After using the quantization, we take the inverse that is defined by the map
We notice that by doing this, the compression becomes lossy during the quantization step as it
creates many zero coefficients as the loss factor becomes larger. As a result, it becomes possible
to spend fewer bits.
3) The functions quantization.m and inverseQuantization.m are using to implement this part of the
lab and are included at the end of the report in section A.2.3 and A.2.4.
4) Below are the six test images reconstructed with a loss factor of 1, 10, and 20 as well as plots of
the reconstruction error, where areas that are more yellow signify a larger error.
6
Figure 1: reconstruction and reconstruction error of barbara.pgm at LF 1
7
Figure 2: reconstruction and reconstruction error of clown.pgm at LF 1
8
Figure 3: reconstruction and reconstruction error of fingerprint.pgm
9
Figure 4: reconstruction and reconstruction error of mandril.pgm with LF 1
10
Figure 5: reconstruction and reconstruction error of roof.pgm with LF1
11
Figure 6: reconstruction and reconstruction error of straw.pgm with LF 1
12
Figure 7: reconstruction and reconstruction error of barbara.pgm with LF 10
13
Figure 8: reconstruction and reconstruction error of clown.pgm with LF 10
14
Figure 9: reconstruction and reconstruction error of fingerprint.pgm with LF 10
15
Figure 10: reconstruction and reconstruction error of mandril.pgm with LF 10
16
Figure 11: reconstruction and reconstruction error of roof.pgm with LF 10
17
Figure 12:reconstruction and reconstruction error of straw.pgm with LF 10
18
Figure 13: reconstruction and reconstruction error of barbara.pgm with LF 20
19
Figure 14: reconstruction and reconstruction error of clown.pgm with LF 20
20
Figure 15: reconstruction and reconstruction error of fingerprint.pgm with LF 20
21
Figure 16: reconstruction and reconstruction error of mandril.pgm with LF 20
22
Figure 17: reconstruction and reconstruction error of roof.pgm with LF 20
23
Figure 18: reconstruction and reconstruction error of straw.pgm with LF 20
24
In figures 1, 7, and 13 the reconstruction is very poor especially around the person as the loss
factor increases. The reconstruction error is quite large for a loss factor of 20, but I cannot really tell the
difference between loss factor of 1 and 10 as I cant really tell different shades of blue very well, but I
decided not to scale these images to tell the reconstruction error better as I dont see the value in it as
the difference between the images doesnt seem to large between loss factors of 1 and 10.
In figures 2, 8, and 14 the story is similar to barbara.pgm as the image has pixels next to each
other are pretty different so the reconstruction is pretty bad for loss factor of 20. For loss factors of 1
and 10 the story is similar to babara.pgm.
In figures 3, 9, and 15 the reconstruction is probably the worst out of all the images as the loss
factor goes up. It is very clear that there are many patches of the fingerprint that are reconstructed
poorly which makes sense as fingerprints have very unique patterns.
In figures 4, 10, and 16 the reconstruction is good around the center of the mandrils face, but
the features around the face are reconstructed poorly as loss factor goes up.
In figures 5, 11, and 17 the reconstruction on the top of the roofs is poor as loss factor goes up.
This makes sense as there are many ridges in the roof that causes the reconstruction to be horrible.
In figures 6, 12, and 18 the reconstruction seems to work reliably and there is only poor
reconstruction for small amounts of pixels rather than a large patch of pixels compared to the other
images when loss factor was high.
5) Below in figure 19 is the PSNR vs. Loss Factor for each of the 6 images. The plots have
logarithmic decay as loss factor increases which makes sense as the reconstruction gets worse
as the loss factor increases and the PSNR has a factor of log 10 in its calculation. It also makes
sense as there is significantly more noise as the reconstruction gets worse which causes the
PSNR to decrease which can be seen in the previous graphs for part 4 as loss factor increases.
25
Figure 19: PSNR vs. Loss Factor
The next part of the lab involves computing the maximum and minimum values of any DCT coefficient or
DC prediction residual can assume. When a quantized coefficient is a zero, the coefficient is not coded.
Instead, we count the number of zeros separating the non-zero coefficient from the previous non-zero
coefficient, we precede it with a code word that encodes the number of zeros between this coefficient
and the previous non-zero coefficients. To achieve this we use arithmetic encoding, and to create the
input to the arithmetic encoder we process the images coefficients into a triplet of numbers: nZeroes,
nBits, and value. These values are based on the DC coefficients of the image and are implemented in the
next assignment.
26
6a) Implementation of the pre-processor run_bits_value can be found at the end of the report in section
A.2.7.
6b) Implementation of the post-processor inverse_run_bits_value is at the end of the report in section
A.2.8.
After implementing the pre-processor, we go to the arithmetic encoding step by feeding in the symb
vector to the arithmetic encoder. This is done in the next assignment.
27
7) The symb vector was produced from run_bits_value and I put it into the irun_bits_value function, but
I did not have time to evaluate the performance of the arithmetic encoder/decoder.
8) This part was skipped as I ran out of time to evaluate the performance of my functions.
3 Conclusion
In this lab, by using the DCT and arithmetic encoding and decoding, we were able to implement
the JPEG algorithm. It could be seen that with greater loss factors that the image reconstruction
becomes worse. By comparing loss factor to the PSNR of the reconstructed images to the original
images as well as the compression ratio to the PSNR, we were able to visually see how well
reconstruction was for certain loss factors.
A MATLAB Code
A.1 Code used to compile report
%This is the main script for lab 5
close all;
clear all;
clc;
pathName = 'C:\Users\house_000\Documents\MATLAB\lab5\images';
fileList = dir(fullfile(pathName,'*.pgm'));
fileName = dir('*.pgm');
numberOfImages = size(fileList,1);
28
for i = 1:numberOfImages
%questions 1-4 parts:
f = imread(fileName(i).name);
figure
imagesc(f, [0 255])
colormap gray
axis square
axis off;
lossFactor = 10;
g = im2double(f);
[coefficientMatrix,temp] = dctmgr(g,lossFactor);
outputImage = idctmgr(coefficientMatrix,temp,lossFactor);
outputImage = im2uint8(outputImage);
figure
imagesc(outputImage, [0 255])
colormap gray
axis square
axis off;
reconstructionError = f - outputImage;
figure
imagesc(outputImage, [0 255])
colormap gray
axis square
axis off;
title(strcat('Reconstructed Image with Loss Factor of: ','
',int2str(lossFactor)));
saveas(gcf,strcat('Image',int2str(i),'
','lossFactor',int2str(lossFactor),'.png'));
figure,imagesc(reconstructionError);
title(strcat('Reconstruction error with Loss Factor of ', '
',int2str(lossFactor)));
saveas(gcf,strcat('Reconstruction error',int2str(i),'
','lossFactor',int2str(lossFactor),'.png'));
end
%question 5
PSNRMatrix = zeros(50,6);
parfor i = 1:numberOfImages
f = imread(fileName(i).name);
%f = imread('mandril.pgm');
g = im2double(f);
for lossFactor = 1:50
lossFactor = 5;
[coefficientMatrix,temp] = dctmgr(g,lossFactor);
symb = run_bits_value(coefficientMatrix);
%minval = min( symb );
%symb = symb - (minval - 1);
%maxval = max( symb );
%histvals = histcounts(symb, maxval);
%histvals = histvals + 1;
%bitstream = arithenco(symb, histvals);
%total_bits = length(bitstream);
%symb_length = length(symb);
%recon_symb = arithdeco(bitstream, histvals, symb_length);
29
%recon_symb = cast(recon_symb, 'double') + double(minval -1);
%irbv_coeffs = irun_bits_value(recon_symb);
outputImage = idctmgr(coefficientMatrix,temp,lossFactor);
outputImage = im2uint8(outputImage);
reconstructionError = f - outputImage;
[peaksnr,snr] = psnr(outputImage, f);
PSNRMatrix(lossFactor,i) = peaksnr;
end
end
figure
lossFactor = 1:50;
fileName = {'barbara', 'clown', 'fingerprint', 'mandril', 'roof', 'straw'}
for i = 1:6
subplot(2,3,i)
plot(lossFactor,PSNRMatrix(:,i));
title(strcat(fileName(i), ' ', 'PSNR vs. Loss Factor'));
end
saveas(gcf, 'LossFactorvsPSNR.png');
A.2 MATLAB functions
A.2.1. dctmgr.m
function [coefficients, temp] = dctmgr(image, lossFactor)
%dctmgr Creates JPEG coefficients for an image
%Takes a luminance (gray-level) image as an input, divives it into
%non overlapping 8x8 blocks, and uses the DCT transform on each block
according
%F(u,v) = 1/4*C_u*C_v sum(f(x,y)*cos(((2x+1)*u*pi)/16) *
%cos(((2y+1)*v*pi)/16), u,v = 0,...,7. For a given block b,
%the coefficients of the column coefficients(:,b) are organized in a
%zig-zag pattern. The zero-frequency coefficients, F(1,1), for each block
%is encoded using differential encoding: coefficients(1,1) = F_1(1,1);
%coefficients(1,2) = F_2(1,1) - F_1(1,1); where F_i(1,1) is the zero
%frequency DCT coefficient of block i
[x,y] = size(image);
index = 1;
coefficients = zeros(64,x*8);
dctMatrix = zeros(x,y);
for x2 = 1:8:x
for y2 = 1:8:y
temp = image(x2:x2+7,y2:y2+7);
firstDct = dct2(temp);
firstDct = quantization(firstDct,lossFactor);
dctMatrix(x2:x2+7,y2:y2+7) = firstDct;
%coefficient matrix:
coefficients(:,index) = zigZag(firstDct);
index = index + 1;
end
end
temp = coefficients(1,:);
for i = 2:y*8
coefficients(1,i) = temp(1,i) - temp(1,i-1);
end
30
A.2.2 idctmgr.m
function [output] = idctmgr(coefficientMatrix,temp,lossFactor)
%(a) It takes a matrix coeff of size 64 x Nblocks, where Nblocks is the
number of 8 x 8 blocks
% inside the image and reconstruct a luminance image.
% (b) For each given block b, the coefficients in the column coeff(:,b) are
used to reconstruct
% the block according to (3).
firstiDct = zeros(8,8);
output = zeros(512,512);
[x,y] = size(coefficientMatrix);
coefficientMatrix(1,:) = temp;
index = 1;
for x2 = 1:8:505
for y2 = 1:8:505
firstiDct = inverseZigZag(coefficientMatrix(:,index));
firstiDct = inverseQuantization(firstiDct,lossFactor);
output(x2:x2+7,y2:y2+7) = idct2(firstiDct);
index = index + 1;
end
end
end
A.2.3 quantization.m
function [output] = quantization(f_uv,lossFactor)
%Implement the quantizer
output = zeros(8,8);
table_Q = [16 11 10 16 24 40 51 61;
12 12 14 19 26 58 60 55;
14 13 16 24 40 57 69 56;
14 17 22 29 51 87 80 62;
18 22 37 56 68 109 103 77;
24 35 55 64 81 104 113 92;
49 64 78 87 103 121 120 101;
72 92 95 98 112 100 103 99];
table_Q = table_Q/255;
for u = 1:8
for v = 1:8
output(u,v) = round((f_uv(u,v)./(lossFactor.*table_Q(u,v))));
end
end
end
A.2.4 inverseQuantization.m
function [output] = inverseQuantization(f_uv,lossFactor)
%Inverse quantization:
output = zeros(8,8);
31
table_Q = [16 11 10 16 24 40 51 61;
12 12 14 19 26 58 60 55;
14 13 16 24 40 57 69 56;
14 17 22 29 51 87 80 62;
18 22 37 56 68 109 103 77;
24 35 55 64 81 104 113 92;
49 64 78 87 103 121 120 101;
72 92 95 98 112 100 103 99];
table_Q = table_Q/255;
for u = 1:8
for v = 1:8
output(u,v) = f_uv(u,v).*lossFactor.*table_Q(u,v);
end
end
end
A.2.5 zigZag.m
function [ outputMatrix ] = zigZag( inputMatrix )
index = reshape(1:numel(inputMatrix),size(inputMatrix));
%flip the matrix left-right. Then use spdiags
%spdiags will get extract the diagional elements and then
%flip the element left-right
index = fliplr(spdiags(fliplr(index)));
%reverse order of odd columns and flip up-down
index(:,1:2:end) = flipud(index(:,1:2:end));
index(index==0) = [];
outputMatrix = inputMatrix(index)';
end
A.2.6 inverseZigZag.m
function [ outputMatrix ] = inverseZigZag(inputMatrix)
%hard coded the inverse ZigZag
% 2. Write the MATLAB function idctmgr that implements the following
functionality:
% (a) It takes a matrix coeff of size 64 x Nblocks, where Nblocks is the
number of 8 x 8 blocks
% inside the image and reconstruct a luminance image.
% (b) For each given block b, the coefficients in the column coeff(:,b) are
used to reconstruct
% the block according to (3).
outputMatrix(1,1) = inputMatrix(1);
outputMatrix(1,2) = inputMatrix(2);
outputMatrix(2,1) = inputMatrix(3);
outputMatrix(3,1) = inputMatrix(4);
outputMatrix(2,2) = inputMatrix(5);
outputMatrix(1,3) = inputMatrix(6);
outputMatrix(1,4) = inputMatrix(7);
32
outputMatrix(2,3) = inputMatrix(8);
outputMatrix(3,2) = inputMatrix(9);
outputMatrix(4,1) = inputMatrix(10);
outputMatrix(5,1) = inputMatrix(11);
outputMatrix(4,2) = inputMatrix(12);
outputMatrix(3,3) = inputMatrix(13);
outputMatrix(2,4) = inputMatrix(14);
outputMatrix(1,5) = inputMatrix(15);
outputMatrix(1,6) = inputMatrix(16);
outputMatrix(2,5) = inputMatrix(17);
outputMatrix(3,4) = inputMatrix(18);
outputMatrix(4,3) = inputMatrix(19);
outputMatrix(5,2) = inputMatrix(20);
outputMatrix(6,1) = inputMatrix(21);
outputMatrix(7,1) = inputMatrix(22);
outputMatrix(6,2) = inputMatrix(23);
outputMatrix(5,3) = inputMatrix(24);
outputMatrix(4,4) = inputMatrix(25);
outputMatrix(3,5) = inputMatrix(26);
outputMatrix(2,6) = inputMatrix(27);
outputMatrix(1,7) = inputMatrix(28);
outputMatrix(1,8) = inputMatrix(29);
outputMatrix(2,7) = inputMatrix(30);
outputMatrix(3,6) = inputMatrix(31);
outputMatrix(4,5) = inputMatrix(32);
outputMatrix(5,4) = inputMatrix(33);
outputMatrix(6,3) = inputMatrix(34);
outputMatrix(7,2) = inputMatrix(35);
outputMatrix(8,1) = inputMatrix(36);
outputMatrix(8,2) = inputMatrix(37);
outputMatrix(7,3) = inputMatrix(38);
outputMatrix(6,4) = inputMatrix(39);
outputMatrix(5,5) = inputMatrix(40);
outputMatrix(4,6) = inputMatrix(41);
outputMatrix(3,7) = inputMatrix(42);
outputMatrix(2,8) = inputMatrix(43);
outputMatrix(3,8) = inputMatrix(44);
outputMatrix(4,7) = inputMatrix(45);
outputMatrix(5,6) = inputMatrix(46);
outputMatrix(6,5) = inputMatrix(47);
outputMatrix(7,4) = inputMatrix(48);
outputMatrix(8,3) = inputMatrix(49);
outputMatrix(8,4) = inputMatrix(50);
outputMatrix(7,5) = inputMatrix(51);
outputMatrix(6,6) = inputMatrix(52);
outputMatrix(5,7) = inputMatrix(53);
outputMatrix(4,8) = inputMatrix(54);
outputMatrix(5,8) = inputMatrix(55);
outputMatrix(6,7) = inputMatrix(56);
outputMatrix(7,6) = inputMatrix(57);
outputMatrix(8,5) = inputMatrix(58);
outputMatrix(8,6) = inputMatrix(59);
outputMatrix(7,7) = inputMatrix(60);
outputMatrix(6,8) = inputMatrix(61);
outputMatrix(7,8) = inputMatrix(62);
outputMatrix(8,7) = inputMatrix(63);
outputMatrix(8,8) = inputMatrix(64);
33
end
A.2.7 run_bits_value.m
function [symb] = run_bits_value(coefficientMatrix)
%implement run_bits_value
nZeros = zeros(64, length(coefficientMatrix));
nBits = zeros(64, length(coefficientMatrix));
zeros_n = 0;
for i = 1:length(coefficientMatrix)
for j = 1:64
value = coefficientMatrix(j,i);
nZeros(j,i) = -1;
if(value == 0)
zeros_n = zeros_n + 1;
end
if j == 1
nBits(j,i) = 12;
else
nBits(j,i) = 11;
end
end
end
symbMatrix = zeros(3*64, length(coefficientMatrix));
mZeros = zeros(64,length(coefficientMatrix));
mBits = zeros(64, length(coefficientMatrix));
mCoefficients = zeros(64,length(coefficientMatrix));
for i = 1:length(coefficientMatrix)
nPresent = 0;
for j = 1:64
if nZeros(j,i) ~= -1
nPresent = nPresent + 1;
mZeros(nPresent,i) = nZeros(j,i);
mBits(nPresent, i) = mBits(j,i);
mCoefficients(nPresent,i) = mCoefficients(j,i);
end
end
end
for i = 1:length(coefficientMatrix)
for j = 1:64
symbMatrix(1+3*(j-1),i) = mZeros(j,i);
symbMatrix(2+3*(j-1),i) = mBits(j,i);
symbMatrix(3+3*(j-1),i) = round(mCoefficients(j,i));
end
end
symbVector = zeros(3*64*length(coefficientMatrix),1);
numI = 1;
for i = 1:length(coefficientMatrix)
for j = 1:64*3 -1
if symbMatrix(j,i) == 0 && symbMatrix(j+1,i) == 0
symbVector(numI) = 0;
symbVector(numI + 1) = 0;
34
numI = numI +2;
break
else
symbVector(numI,1) = symbMatrix(j,i);
numI = numI + 1;
end
end
end
numI = 0;
for i = 1:length(symbVector) - 3
if symbVector(i) == 0 && symbVector(i+1) == 0 && symbVector(i+2) ==
0 && symbVector(i+3) == 0
break;
else
numI = numI +1;
end
end
symb = symbVector(1:numI);
end
A.2.8 inverse_run_bits_value.m
function [ coefficients ] = inverse_run_bits_value( symb )
%%(6)
blockN = 1;
start = zeros(4096,1);
start(1) = 1; %first block starts at 1
ends(4096) = length(symb); %ast block ends at the last symb element
for i = 3:length(symb)-1 %starting index of every block
if symb(i-2) == 0 && symb(i-1) == 0 && symb(i) == 0 %find EOB when three
consecutive zeros
start(blockN+1) = i;
blockN = blockN +1;
end
end
matComp = zeros(64,4096); %compressed matrix, doing decompression
for i = 1:4095
matComp(1:length(symb(start(i):start(i+1)-1)),i) =
symb(start(i):start(i+1)-1);
end
%handle last entry separately because of bound conditions
matComp(1:length(symb(start(4096):length(symb))),4096) =
symb(start(4096):length(symb));
matUncomp = zeros(64,4096); %uncompressed matrix
for i = 1:4096
ind = 0; %keep track of where the value is
for j = 1:3:61 %go in jumps of 3
zToSkip = matComp(j,i); %get the number of zeros that need to be
placed in between
nBits = matComp(j+1,i);
val = matComp(j+2,i); %get the value from the triplet
35
ind = ind+zToSkip+1; %new index to place val is the old
index+numZeros+1
matUncomp(ind,i) = val;
end
end
coefficients = matUncomp;